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I .  INTRODUCTION 

-  The  study  of  technical  solutions  to  the  problem  of  managing  fluid 
waste  discharges  in  an  estuary  involves  complicated  relationships,  such 
as  those  between  the  waste  load,  the  location  of  discharges,  the  degree 
of  treatment,  the  geometry  of  the  estuary,  the  flow  in  the  estuary,  and 
the  temperature.  Models  can  be  built  to  predict  the  probable  conse¬ 
quences  of  each  of  the  alternative  solutions  chosen.  Using  the  pre¬ 
dictions  obtained  from  these  models  and  whatever  other  information  or 
insight  is  relevant,  alternatives  can  be  compared  and  conclusions  drawn 
concerning  the  most  desirable  course  of  action.  Such  a  model  for  well- 
mbted  estuaries  and  coastal  seas  is  presented  Kara  and  i»  described. ip 
more  detail  in  Reference  1.  The  model  is  being  used  and  evaluated  in 

a  study  of  the  extent  and  control  of  pollution  in  Jamaica  Bay,  Npw 

(2) 

York.  Because  the  fluid  waste  discharged  into  the  estuary  is  trans¬ 
ported  in  it,  a  model  of  the  fluid  must  be  made  in  order  to  establish 
the  relationship  between  the  transport  of  the  pollution  and  the  trans¬ 
porting  medium.  Also,  it  can  be  visualized  that  constituents  are  trans¬ 
formed  during  the  transport  by  chemical  or  biochemical  action,  and  for 
this  process  a  reaction  model  must  be  designed. 

The  flow  in  the  estuary  is  governed  by  the  shape  of  the  estuary 
and  its  boundaries  at  the  sea  and  the  inland  ends,  as  well  as  by  the 
tides  and  the  river  discharges.  The  latitude  of  the  area  under  inves¬ 
tigation  determines  the  effect  of  the  earth's  rotation  on  the  flow, 
which  is  of  particular  importance  in  large  estuaries  and  coastal  seas 
at  higher  latitudes. 

In  the  flow  model,  velocities  of  the  water  in  the  region  are  com¬ 
puted  together  with  the  water  levels.  Because  the  tide  stage  influences 
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the  location  of  the  water's  edge  on  the  tidal  flats,  which  affects  sur¬ 
face  area  and  cross  sections  to  a  considerable  extent,  such  functional 
relationships  must  also  be  included  in  the  model. 

The  fluid  waste  loads  in  the  system  influence  the  pollution  con¬ 
stituent  transport,  while  air  temperature,  solar  radiation,  and  ex¬ 
change  of  gases  at  the  water  surface  Influence  the  reactions. 

In  engineering  and  in  systems  analysis  relationships  are  custom¬ 
arily  displayed  in  graphical  form.  Because  of  the  complexity  of  *e 
results  obtained  from  the  simulations,  the  graphical  results  produced 
by  the  model  permit  an  insight  into  the  movements  of  the  pollutants  in 
the  estuary  in  relation  to  the  flow.  The  procedures  used  to  derive 
these  graphs  are  an  integral  part  of  the  simulation  model. 


II .  MOMENTUM  AND  MASS  TRANSPORT 


The  fluid  flow  in  an  estuary  is  three-dimensional  and  time-depen¬ 
dent.  Because  the  flow  is  mainly  horizontal,  the  classic  approach  is 
to  assume  that  the  pressures  are  hydrostatic  and  that  only  shear  stresses 
from  horizontal  velocity  components  are  important. 

In  this  classical  approach,  it  is  further  assumed  that  gross  mo¬ 
tions  only  ere  considered  and  the  effects  of  small-scale  velocity 
fluctuation  are  aggregated  into  shear  s tress  terms. 

At  present  (1969) ,  computational  techniques  are  Inadequate  to 
deal  with  three-dimensional  computations  of  fluid  flow,  and  problems 
with  the  complicated  boundaries  of  estuaries  are  beyond  the  capabil¬ 
ities  of  present  computers.  In  the  approach  presented  here,  vertical 
integration  of  the  equation  of  motion  and  continuity  is  used  to  reduce 

the  problem  to  a  two-dimensional  one. 

(1  3) 

Such  integration  yields'  * 
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where  U,V  are  the  over-che-vertical  averaged  velocities  along  the  east¬ 
ward  and  northward  directed  axes,  respectively,  p  is  the  density  of 
water,  g  the  gravity  acceleration,  f  the  Coriolis  parameter,  t  the  time, 

H  the  temporal  depth,  i;  the  waterlevel  elevation  relative  to  the  ref- 

x  y 

erence  plane,  C  the  Chezy  term  and  t  and  t7  the  windstress  components 
on  a  finite  surface  element. 

Vertical  integration  of  the  mass  balance  equations  of  dilute  pol¬ 
lution  constituents  in  well-mixed  estuaries  results  in  the  expressions 
for  the  transport  model 

3 (HP)  9 (HUP)  3(HVP)  3(HDx  3^)  3(HDy  3y)  _  _  n 

3t  3x  3y  3x  3y  U  w 

where  P  is  the  vector  of  mass  concentration  of  constituents  averaged 

over  the  temperal  depth,  D  ,D  are  the  dispersion  coefficient  and  R  is 

x  y 

a  vector.  The  elements  of  this  vector  include  the  local  addition  of 
the  constituents  in  addition  to  the  rate  of  production  of  the  constit¬ 
uent  in  a  water  column  w‘-h  unit  dimension. 

For  any  point  in  the  model,  the  following  functional  relation 
is  used: 


f1(U,C,H)  +  Dw 

(5) 

f2(V,C,H)  +  Dw 

(6) 

where  Dw  is  a  diffusion  coefficient  dependent  on  wave  and  wind  condi¬ 
tions  and  the  lateral  diffusion. 


III.  THE  COMPUTATIONAL  MODEL 

The  discrete  values  of  the  variables  are  described  on  a  grid.  A 
space-staggered  scheme  is  used,  in  which  water  levels  and  velocities 
are  described  at  different  grid  points  (Fig.  1).  The  mass  densities 
are  described  at  the  same  location  as  the  water  levels.  This  scheas 
has  the  advantage  that  in  the  formula  operated  upon  in  time  there  exists 
a  centrally  located  spatial  derivative  for  the  linear  tern.  The  water 
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+  Water  level  (.f)  and  mass  density  p 

o  Depth  (h) 

-  U  velocity  (u) 

I  V  velocity  (v) 

x 

Fig.  1— Space-staggered  scheme 

level  r,  and  the  mass  densities  are  described  at  integer  values  of 
j  and  k,  the  velocities  u  are  described  at  half-integer  values  of  j 
and  integer  values  of  R,  and  the  velocity  v  is  described  at  Integer 
values  of  j  and  half-integer  values  of  k.  The  water  depth  is  described 
at  the  half-integer  values  of  j  and  k. 

Fractional  and  whole  time  steps  are  used:  The  fractional  steps 
are  used  for  the  equations  of  continuity  for  the  fluid  mass  and  the 
mass  densities,  and  the  whole  steps  for  the  velocity  componeats  u  and 
v.  The  velocity  components  are  staggered  in  time. 

The  following  notation  is  used  for  the  approximation  of  the  dif¬ 
ferential  equations  by  a  system  of  difference  equations: 

F  a  F(JAx,  kAy,  nAt) 

where  (x,y,t)  •  (JAx,  kAy,  nAt)  and  j,k,t  ■  0,  ±$,  ±1,  ±3/2,  ... 

The  following  symbolic  expressions  are  used  for  averages  and  dif¬ 
ferences.  They  are  shorn  only  for  x,  but  are  also  used  for  y  and  t: 

F*  •  i|F{(J+£)  Ax, kAy, nAt*  ♦  F[(J-$)  Ax,kAy,nAt]j  (7) 

{xF  *  Ax, kAy, nAt]  -  F[(j-$)  Ax,kAy,nAt]| 


(8) 
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F  »  \  F[ (J+J)  Ax,  (k+i)  Ay.nAt]  +  F[(j+£)  Ax,  (k-£  Ay,nAt] 

+  F[(j-£)  Ax,  (k4-£)  Ay,nAt]  +  F((j-£)  Ax,  (k-§)  Ay.nAt]  (9) 

Special  notations  are  used  only  for  the  Indication  of  shifted  time 
levels . 

6+ft  F  *  ft{F[JAx,kAy,  (n+£)  At]  -  F(jAx,kAy,nAt)j  (10) 

F+  -  F[jAx,kAy,  (n+£)  At]  (11) 

F_  -  FfjAx.kAy,  (n-i)  At]  (12) 


The  sets  of  finite-difference  equations  will  now  be  presented  in 
two  sequential  groups  as  they  are  used,  and  their  use  will  be  discussed 
subsequently.  For  the  time  level  n,  we  use 


5  u  -  fv  +  u  6  uM  +  v  6  uy  +  g  -  g 

t  +  x  -  y  •  x*  • 


?  fu  >2  +  ci)2]* 


(i?  +  7>)(c*)2 

at  J+$,k,n  (13) 


at  J ,k,n 


at  J ,k,n 

[«« *  o]  *  *,  [(p  *  ?)  vS]  * #  *  ?)  vf] 


♦  R  •  0 
at  J,k,n 


and  the  finite  difference  approximations  at  time  level  n4$  ere 
*  v . ^  t  ?[<5>2  ♦  (v.)2]* 

*«’*fc*‘  ij?*.  V  •«  (r-4?)^)2  '° 

et  J.k^.tH*  (16) 
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S4*t5  +  6x 


(?  ♦  ?)  »]  +  «,[(>?  *  F)  vj  -  o 


at  j ,k,n+£ 

6Ht[p(F +  0] + »,[(s' +  ^ u?] + *,[(* +  F)  v^l 


+  -A  d  6  p 

-  <5  (ii*  +  D  6  P. 

L\  /  x  X 

y  [V  +/  y  y  +_ 

(17) 


(18) 


The  sets  of  equations  can  be  solved  very  rapidly  as  described  in 
Ref.  1.  The  method  consists  of  solving  a  set  of  linear  equations,  each 
having  three  unknown  values. 


IV.  SIMULATIONS 

Jamaica  Bay  is  the  area  for  which  the  model  is  applied.  It  has 
a  very  complicated  bathymetry;  many  natural  and  dredged  channels  are 
separated  by  tidal  flats.  Thus,  a  finite-difference  approximation  of 
the  bay  bottom  with  a  small  grid  size  is  required.  The  number  of  grid 
points  upon  which  the  computations  are  made,  however,  is  limited  by 
the  size  of  the  memory  of  presently  available  computers. 

A  grid  size  of  500  ft  was  selected.  This  grid  size  permitted  sev¬ 
eral  grid  points  to  be  present  across  the  main  channels  in  the  bay  and 
left  space  in  the  computer  memory  for  program  growth  during  development. 
The  computational  arrays  were  78  *  61  points.  Two  arrays  each  are  used 
for  the  two  velocity  components  and  three  arrays  for  the  water  lovels, 
since  three  time  levels  are  involved.  One  array  is  used  for  depth  data 
and  one  for  the  Chezy  coefficients.  If  three  pollution  substances  are 
computed,  six  arrays  are  necessary,  i.e.,  one  for  each  of  the  subse¬ 
quent  time  levels,  thus,  in  total,  15  arrays  are  to  be  kept  in  memory, 
which  requires  approximately  75,000  storage  locations. 

Since  the  problem  is  time-dependent,  large  volumes  of  input  data 
are  generally  required,  these  input  data  consist  of  the  time  histories 
of  the  following: 

1.  Water  levels  end  mass  densities  at  the  mouth  of  the  bay. 

2.  Discharges  at  the  location  o*  the  outfalls  of  treatment  plants 
S' d  the  combined  sever  over i lows. 
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3.  Mass  densities  associated  with  these  discharges. 

4.  Wind  speed  and  direction. 

The  environment  simulation  program  described  in  tula  paper  gener¬ 
ates  large  volumes  of  data  in  numerical  form.  Printed  numerical  data 
become  unmanageable,  and  extracting  important  data  from  a  large  volume 
of  material  is  a  considerable  task.  Thus,  programs  were  written  for 
machine-made  drawings  to  represent  conditions  at  selected  times. 

f  3) 

Using  the  Integrated  Graphics  System,  which  supports  an  S-C  4060, 
gra.hic  outputs  were  made  of  the  type  shown  in  Fig.  2. 

This  graphic  output  displays  the  outline  of  the  bay  (approximate 
location  of  the  high  water  line)  and  the  grid  points  upon  which  the 
water  levels  and  mass  densities  are  computed.  The  total  number  of  points 
varies  from  output,  to  output  as  the  locations  of  the  water's  edge  in  the 
marshy  areas  determine  which  point  is  included.  The  blank  areas  with 
no  dots  (.)  are  thus  indications  for  land. 

The  U  velocity  component  is  computed  at  a  location  between  the 
water  level  points  in  the  x-direction,  and  the  V  velocity  component  is 
computed  at  a  location  between  water  level  points  In  the  y-dlrectlon 
(Fig.  1).  In  order  to  display  the  magnitude  and  direction  of  the  ve¬ 
locity  in  the  field  of  computation,  velocity  vectors  were  computed  at 
the  locations  of  the  water  levels  from  the  average  values  of  the  adja¬ 
cent  components.  To  enhance  the  representation,  vectors  are  plotted 
only  at  every  second  point  In  the  x-  and  y-directlnns. 

Variables  such  as  mass  densities  of  constituents  (or  the  water  lev¬ 
els)  are  plotted  by  means  of  isocontour  lines.  Different  line  thick¬ 
ness  or  dashing  Is  available  for  such  presentation.  The  subroutines 
to  draw  such  contours  are  very  extensive,  since  a  search  must  be  made 
at  each  location  if  the  contour  inder  consideration  passes  through  the 
adjacent  grids.  In  addition,  proper  terminations  of  such  contours  had 
to  be  made  If  the  adjacent  grid  points  had  become  land  and  were  thus 
not  computed. 

The  value  of  each  lsocontour  line  can  be  found  from  extrapolation 
of  numerical  values  which  are  also  plotted  In  the  field  at  selected  lo¬ 
cations. 
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Figure  2  presents  the  pattern  of  the  distribution  of  coliform  bac' 
teria  from  a  source  in  the  northwest  section  of  the  bay  8  hours  after 
the  start  of  the  computation. 

The  model  is  now  being  verified  for  its  hydraulic  and  water  al- 
ity  properties.  The  verification  is  based  upon  an  extensive  field  in¬ 
vestigation  and  analysis. 
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